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Abstract 

We  present  a  generalized  polynomial  chaos  algorithm  to  model  the  input  uncertainty  and  its  propagation  in 
flow-structure  interactions.  The  stochastic  input  is  represented  spectrally  by  employing  orthogonal  polynomial 
functionals  from  the  Askey  scheme  as  the  trial  basis  in  the  random  space.  A  standard  Galerkin  projection  is 
applied  in  the  random  dimension  to  obtain  the  equations  in  the  weak  form.  The  resulting  system  of  deterministic 
equations  is  then  solved  with  standard  methods  to  obtain  the  solution  for  each  random  mode.  This  approach  is 
a  generalization  of  the  original  polynomial  chaos  expansion,  which  was  first  introduced  by  N.  Wiener  (1938)  and 
employs  the  Hermite  polynomials  (a  subset  of  the  Askey  scheme)  as  the  basis  in  random  space.  The  algorithm  is 
first  applied  to  second-order  oscillators  to  demonstrate  convergence,  and  subsequently  is  coupled  to  incompressible 
Navier-Stokes  equations.  Error  bars  are  obtained,  similar  to  laboratory  experiments,  for  the  pressure  distribution 
on  the  surface  of  a  cylinder  subject  to  vortex-induced  vibrations. 
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1  Introduction 


In  the  last  decade  there  has  been  substantial  progress  in  simulations  of  flow-structure  interactions  involving  the  full 
Navier-Stokes  equations,  e.g.  [1,  2].  While  such  simulations  are  useful  in  complementing  experimental  studies  in 
the  low  Reynolds  number  range,  they  are  based  on  ideal  boundary  conditions  and  precisely  defined  properties  of  the 
structure.  In  practice,  such  flow  conditions  and  properties  can  only  be  defined  approximately.  As  an  example,  the 
internal  structural  damping  for  the  structure  is  typically  taken  as  1  —  3%  of  the  critical  damping  since  it  cannot  be 
quantified  by  direct  measurements.  It  is,  therefore,  of  great  interest  to  formally  model  such  uncertainty  of  stochastic 
inputs,  and  to  formulate  algorithms  that  reflect  accurately  the  propagation  of  this  uncertainty  [3] . 

To  this  end,  the  Monte  Carlo  approach  can  be  employed  but  it  is  computationally  expensive  and  is  only  used  as 
the  last  resort.  The  sensitivity  method  is  a  more  economical  approach,  based  on  the  moments  of  samples,  but  it  is 
less  robust  and  depends  strongly  on  the  modeling  assumptions  [4] .  One  popular  technique  is  the  perturbation  method 
where  all  the  stochastic  quantities  are  expanded  around  their  mean  via  Taylor  series.  This  approach,  however,  is 
limited  to  small  perturbations  and  does  not  readily  provide  information  on  high-order  statistics  of  the  response.  The 
resulting  system  of  equations  becomes  extremely  complicated  beyond  second-order  expansion.  Another  approach 
is  based  on  expanding  the  inverse  of  the  stochastic  operator  in  a  Neumann  series,  but  this  too  is  limited  to  small 
fluctuations,  and  even  combinations  with  the  Monte  Carlo  method  seem  to  result  in  computationally  prohibitive 
algorithms  for  complex  systems  [5] . 

A  more  effective  approach  pioneered  by  Ghanem  &  Spanos  [6]  in  the  context  of  finite  elements  for  solid  mechanics 
is  based  on  a  spectral  representation  of  the  uncertainty.  This  allows  high-order  representation,  not  just  first-order  as 
in  most  perturbation-based  methods,  at  high  computational  efficiency.  It  is  based  on  the  original  theory  of  Wiener 
(1938)  on  homogeneous  chaos  [7,  8].  This  approach  was  employed  in  turbulence  in  the  1960s  [9,  10,  11].  However, 
it  was  realized  that  the  chaos  expansion  converges  slowly  for  turbulent  fields  [12,  13,  14],  so  the  polynomial  chaos 
approach  did  not  receive  any  attention  for  a  long  time. 

In  more  recent  work  [15,  16]  the  polynomial  chaos  concept  was  extended  to  represent  many  different  distribution 
functions.  This  generalized  polynomial  chaos  approach,  also  referred  as  the  Askey-chaos,  employs  the  orthogonal 
polynomials  from  the  Askey  scheme  [17]  as  the  trial  basis  in  the  random  space.  The  original  polynomial  chaos  can 
be  considered  as  a  subset  of  the  generalized  polynomial  chaos,  as  it  employs  Hermite  polynomials,  a  subset  of  the 
Askey  scheme,  as  the  trial  basis.  In  [15],  the  framework  of  Askey-chaos  was  proposed  and  convergence  properties  of 
different  random  bases  were  examined.  In  [16]  the  Askey-chaos  was  applied  to  model  uncertainty  in  incompressible 
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Navier-Stokes  equations.  Various  tests  were  conducted  to  demonstrate  the  convergence  of  the  chaos  expansion  in 
prototype  flows. 

For  flow-structure  interactions  the  interest  on  stochastic  modeling  so  far  has  primarily  been  on  the  dynamics  of 
lumped  systems,  i.e. ,  single-  or  two-degree-of- freedom  second-order  oscillators  [18,  19].  The  effect  of  the  flow  has 
been  modeled  via  an  interaction  (source)  term  as  either  white  noise  or  as  a  Gaussian  distribution  if  the  loading  is 
caused  by  wind  [20,  21,  22,  23].  However,  non-Gaussian  distribution  behavior  for  the  response  has  been  documented 
with  the  excess  index  well  above  or  below  zero  (sharp  or  flat  intermittency)  [18].  For  example,  even  for  a  velocity 
field  following  a  Gaussian  distribution,  which  is  a  reasonable  assumption  for  maritime  winds  [21],  the  corresponding 
force  given  by  the  Morison  formula 

Fv{t)=l-PDCdV{t)\V{t)\ 


does  not  follow  a  Gaussian  distribution.  This  is  because  the  above  formula  defines  a  nonlinear  (memoryless)  trans¬ 
formation  [20] ,  and  its  first-density  function  is  given  by 


1  sign (v)^/\v\-  mv  2 


(JV 


n 


where  my  and  ay  are  the  mean  value  and  variance  of  the  the  Gaussian  distribution  for  the  velocity  V(t). 


In  this  paper  we  apply  the  chaos  expansions  to  coupled  Navier-Stokes/structure  equations.  We  first  demonstrate 
the  convergence  of  chaos  expansions  by  solving  a  second-order  ordinary  differential  equation.  We  then  present  the 
stochastic  modeling  of  the  fully  coupled  flow-structure  interaction  problem  for  vortex-induced  vibrations  in  flow  past 
a  cylinder.  The  algorithms  developed  here  are  general  and  can  be  applied  to  any  type  of  distributions  although  our 
applications  are  concentrated  on  Gaussian  type  random  inputs. 

In  the  next  section  we  review  the  theory  of  the  generalized  polynomial  chaos.  In  section  3  we  apply  it  to  second- 
order  oscillators,  and  in  section  4  we  present  its  application  to  Navier-Stokes  equations.  In  section  5  we  present  the 
computational  results  of  stochastic  flow-structure  interactions,  and  we  conclude  with  a  brief  discussion  in  section  6. 

2  The  Generalized  Polynomial  Chaos 

In  this  section  we  introduce  the  generalized  polynomial  chaos  expansion  along  with  the  Karhunen-Loeve  (KL)  ex¬ 
pansion,  another  classical  technique  for  representing  random  processes.  The  KL  expansion  can  be  used  in  some  cases 
to  represent  efficiently  the  known  stochastic  fields,  i.e.,  the  stochastic  inputs. 
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2.1  The  Askey  Scheme 


The  Askey  scheme,  which  is  represented  as  a  tree  structure  in  figure  1  (following  [24]),  classifies  the  hypergeometric 
orthogonal  polynomials  and  indicates  the  limit  relations  between  them.  The  ‘tree’  starts  with  the  Wilson  polynomials 
and  the  Racah  polynomials  on  the  top.  The  Wilson  polynomials  are  continuous  while  the  Racah  polynomials  are 
discrete.  The  lines  connecting  different  polynomials  denote  the  limit  transition  relationships  between  them;  this 
implies  that  the  polynomials  at  the  lower  end  of  the  lines  can  be  obtained  by  taking  the  limit  of  one  of  the  parameters 
from  their  counterparts  on  the  upper  end.  For  example,  the  limit  relation  between  Jacobi  polynomials  (x)  and 

Hermite  polynomials  Hn(x)  is 


lim  a-*"p(«>“)  (4=1  = 


Hn(x) 

2 nn\  ’ 


and  between  Meixner  polynomials  Mn(x:  (3,  c)  and  Charlier  polynomials  Cn(x;a)  is 


lim  Mn  (  x;  (3,  Q  )  =  Cn(x ;  a). 

0—>oo  \  a  +  p ) 

For  a  detailed  account  of  definitions  and  properties  of  hypergeometric  polynomials,  see  [17];  for  the  limit  relations  of 
Askey  scheme,  see  [25]  and  [24]. 


Figure  1:  The  Askey  scheme  of  orthogonal  polynomials 

The  orthogonal  polynomials  associated  with  the  generalized  polynomial  chaos,  include:  Hermite,  Laguerre,  Jacobi, 
Charlier,  Meixner,  Krawtchouk  and  Hahn  polynomials. 
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2.2  The  Generalized  Polynomial  Chaos:  Askey-Chaos 

The  originial  polynomial  chaos  [7,  8]  employs  the  Hermite  polynomials  in  the  random  space  as  the  trial  basis  to 
expand  the  stochastic  processes.  Cameron  &  Martin  proved  that  such  expansion  converges  to  any  second-order 
processes  in  the  L2  sense  [26].  It  can  be  seen  from  figure  1  that  Hermite  polynomial  is  a  subset  of  the  Askey  scheme. 
The  generalized  polynomial  chaos,  or  the  Askey-Chaos,  was  proposed  in  [15,  16]  and  employs  more  polynomials  from 
the  Askey  scheme.  Convergence  to  second-order  stochastic  processes  can  be  readily  obtained  as  a  generalization  of 
Cameron-Martin  theorem  [26]. 

A  general  second-order  random  process  A(0),  viewed  as  a  function  of  9  €  (0,1),  i.e.  the  random  event,  can  be 
represented  in  the  form 


X(0)  =  a0I0 

OO 

+  ]T  ^Jr^)) 

*1=1 
00  i\ 

+  X!  Ciii2-f2(61(6,),^2(6')) 

il  —  1  22  —  1 
OO  %\  22 

+  EEE 

*1=1  *2=1  *3  =  1 

+  • • •  ,  (1) 

where  ■■■)£*„)  denotes  the  Askey-chaos  of  order  n  in  terms  of  the  multi-dimensional  random  variables  £  = 

(£n ,  In  the  original  polynomial  chaos,  {/„}  are  Hermite  polynomials  and  £  are  Gaussian  random  variables. 

In  the  Askey-chaos  expansion,  the  polynomials  In  are  not  restricted  to  Hermite  polynomials  and  £  not  Gaussian 
variables.  The  corresponding  type  of  polynomials  and  their  associated  random  variables  are  listed  in  table  1. 


Random  variables  £ 

Orthogonal  polynomials  { /„  } 

Support 

Continuous 

Gaussian 

Hermite 

Gamma 

Laguerre 

Beta 

Jacobi 

[a,b\ 

Uniform 

Legendre 

[a,b] 

Discrete 

Poisson 

Charlier 

{0,1,2,...} 

Binomial 

Krawtchouk 

{0,1, 

Negative  Binomial 

Meixner 

{0,1,2,...} 

Hypergeometric 

Hahn 

{0,1, ...,7V} 

Table  1:  Correspondence  of  the  type  polynomials  and  random  variables  for  different  Askey-chaos  (TV  >  0  is  a  finite 
integer). 


5 


For  notational  convenience,  we  rewrite  equation  (1)  as 

OO 

*(*)  =  £W*),  (2) 

j= o 

where  there  is  a  one-to-one  correspondence  between  the  functions  Ini^h ,  •  •  • ,  £in)  and  <Fj(£),  and  their  coefficients 
Cj  and  Since  each  type  of  polynomials  from  the  Askey  scheme  form  a  complete  basis  in  the  Hilbert  space 

determined  by  their  corresponding  support,  we  can  expect  each  type  of  Askey-chaos  to  converge  to  any  L2  functional 
in  the  L2  sense  in  the  corresponding  Hilbert  functional  space  as  a  generalized  result  of  Cameron-Martin  theorem 
([26]  and  [27]).  The  orthogonality  relation  of  the  Askey-Chaos  polynomial  chaos  takes  the  form 

<  >=<  $2  >  Sijj  (3) 


where  Si7  is  the  Kronecker  delta  and  <  •,  •  >  denotes  the  ensemble  average  which  is  the  inner  product  in  the  Hilbert 
space  of  the  variables  £ 

<  mate)  >=  j  mg(®w(Z)dZ,  (4) 

or 

</(€)$(€)>=  (5) 

€ 

in  the  discrete  case.  Here  W(£)  is  the  weighting  function  corresponding  to  the  Askey  polynomials  chaos  basis 
{4>i}.  Each  type  of  orthogonal  polynomials  from  the  Askey-chaos  has  weighting  functions  of  the  same  form  as  the 
probability  function  of  its  associated  random  variables  £,  as  shown  in  table  1. 

For  example,  as  a  subset  of  the  Askey-chaos,  the  original  polynomial  chaos,  also  will  be  termed  the  Hermite-chaos, 
employs  the  Hermite  polynomials  defined  as 


4(&,.  ■•,&„)  =ei*  €(-l) 


dn 


d&i  •  •  •  d£in 


.-*eTe 


(6) 


where  £  =  are  multi-dimensional  independent  Gaussian  random  variables  with  zero  mean  and  unit 

variance.  The  weight  function  in  the  orthogonality  relation  (4)  is 


TT(£) 


V(2ir)r 


(7) 


where  n  is  the  dimension  of  It  can  seen  that  this  is  the  same  as  the  probability  density  function  (PDF)  and  n- 
dimensional  Gaussian  random  variables.  For  example,  the  one-dimensional  Hermite  polynomials  are: 


*0  =  1,  *!  =  ■£,  *2=£2-l,  *3  =  £3  -  3£, 


(8) 
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2.3  The  Karhunen-Loeve  Expansion 


The  Karhunen-Loeve  (KL)  expansion  [28]  is  another  way  of  representing  a  random  process.  It  is  based  on  the  spectral 
expansion  of  the  covariance  function  of  the  process.  Let  us  denote  the  process  by  /i(x,  9)  and  its  covariance  function 
by  where  x  and  y  are  the  spatial  or  temporal  coordinates.  By  definition,  the  covariance  function  is  real, 

symmetric,  and  positive  definite.  All  eigenfunctions  are  mutually  orthogonal  and  form  a  complete  set  spanning  the 
function  space  to  which  /i(x,  9)  belongs.  The  KL  expansion  then  takes  the  following  form: 

OO 

Mx,  0)  =  h(x)  +  ^2  \AI<Mx)&(#)>  (9) 

i=l 

where  /i(x)  denotes  the  mean  of  the  random  process,  and  £,i{9)  forms  a  set  of  independent  random  variables.  Also, 
0i(x)  and  \  are  the  eigenfunctions  and  eigenvalues  of  the  covariance  function,  respectively,  i.e. , 

J  Rhh{x ,  y)<My)dy  =  A ^(x).  (10) 

Among  many  possible  decompositions  of  a  random  process,  the  KL  expansion  is  optimal  in  the  sense  that  the 
mean-square  error  of  the  finite  representation  of  the  process  is  minimized.  Its  use,  however,  is  limited  as  the  covariance 
function  of  the  solution  process  is  often  not  known  a  priori.  Nevertheless,  the  KL  expansion  provides  an  effective 
means  of  representing  the  input  random  processes  when  the  covariance  structure  is  known. 


3  Second-order  Random  Oscillator 

3.1  Governing  Equations 


We  consider  the  second-order  linear  ordinary  differential  equation  (ODE)  system  with  both  external  and  parametric 
random  excitations. 


dx 

dt 


=  y , 


(ii) 


[  ^ +c(9)y  +  k(9)x  = 

where  the  parameters  and  forcing  are  functions  of  random  event  8.  We  assume 

c=c  +  ctc£  i,  k  =  k  +  Gkt.2,  /(f)  =  Fcos(wi)  =  (/+ 0763)  cos(wt),  (12) 

where  (c,  ac),  ( k,ak )  &nd  (f^f)  are  mean  and  standard  deviation  of  c,  k  and  F,  respectively.  The  random 
variables  £i,  £2  and  £3  are  assumed  to  be  independent  standard  Gaussian  random  variables. 


3.2  Chaos  Expansions 


By  applying  the  generalized  polynomial  chaos  expansion,  we  expand  the  solutions  as 

p  p 

x(t)  =  y(t )  = 

i— 0  i— 0 


(13) 
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where  we  have  replaced  the  infinite  summation  in  infinite  dimension  of  C  in  equation  (2)  by  a  truncated  finite- 
term  summation  in  finite  dimensional  space  of  C-  In  this  case,  £  =  (£1,^2,  £3)  is  a  three-dimensional  Gaussian 
random  vector  according  to  the  random  inputs.  This  results  in  a  three-dimensional  Hermite- chaos  expansion.  The 
most  important  aspect  of  the  above  expansion  is  that  the  random  processes  have  been  decomposed  into  a  set  of 
deterministic  functions  in  the  spatial-temporal  variables  multiplied  by  the  random  basis  polynomials  which  are 
independent  of  these  variables: 


Ef=0  =  £f=0 


Efc=o  dt  +  Ei=0  E;=o  +  Ei=0  E,'=0  —  Efc=o  fk{t)&k> 


(14) 


where  Cj,  ki  and  ft  are  the  chaos  expansion,  similar  to  equation  (13),  of  c,  k  and  /,  respectively.  A  Galerkin  projection 
of  the  above  equation  onto  each  polynomial  basis  {4b}  is  then  conducted  in  order  to  ensure  the  error  is  orthogonal 
to  the  functional  space  spanned  by  the  finite-dimensional  basis  {4>;}.  By  projecting  with  4>fc  for  each  k  =  {0, . . . ,  P} 
and  employing  the  orthogonality  relation  (3),  we  obtain  for  each  k  =  0, . . .  P, 


dt 


(15) 


dt  <$£>  Ei=o  Ej=o(c*2/i  +  kiXj)eijk  —  . fk(t ), 
where  e^-fc  =<  4>i4>:)4>fc  >.  Together  with  <  >,  the  coefficients  can  be  evaluated  analytically  from  the 

definition  of  4q.  Equation  (15)  is  a  set  of  (P  +  1)  coupled  ODEs.  The  total  number  of  equation  is  determined  by  the 
dimensionality  of  the  chaos  expansion  (n),  in  this  case  (n  =  3),  and  the  highest  order  (p)  of  polynomials  4>  [6]: 


P  1  s-l 


p=E^ri(n+r)- 


S! 

s=l  r— 0 


(16) 


3.3  Numerical  Results 


The  above  set  of  equations  can  be  integrated  by  any  conventional  method,  e.g.,  Runge-Kutta.  Here  we  employ  the 
Newmark  scheme  which  is  second-order  accurate  in  time.  We  define  two  error  measures  for  the  mean  and  variance 
of  the  solution 


a(T)  = 


X(T)  -  Zexact  (T) 


d'exrict.  (t  ) 


£va,r(T)  = 


\T )  -  , 


tCH 


t(T) 


(17) 


where  x(t)  =  E[x(t)]  is  the  mean  value  of  x(t)  and  er2(f)  =  E  ^(a’(f)  —  ir(f))2J  is  the  variance.  Integration  is  performed 
up  to  T  =  100  (nondimensional  time  units)  when  the  solution  reaches  an  asymptotic  periodic  state.  The  computation 
parameters  are  set  as:  (c,  trc)  =  (0.1,0.01),  (fc,crfc)  =  (1.05,0.105)  and  (/,  07)  =  (0.1,0.01),  with  frequency  to  =  1.05 
and  zero  initial  conditions.  The  exact  stochastic  solution  is  obtained  from  the  exact  deterministic  solution  and 


the  known  probability  distribution  functions  of  the  random  inputs.  We  further  need  to  integrate  the  solution  over 


the  support  defined  by  the  Gaussian  distribution  to  obtain  the  exact  mean  and  variance  of  the  solution.  These 
integrations  are  performed  numerically  using  a  Gauss-Hermite  quadrature;  a  quadrature  with  30  points  provides 
high  accuracy. 

In  figure  2  (left)  we  plot  the  development  of  mean  (zero  mode)  as  well  as  the  first  three  random  modes,  i.e.  the 
modes  contribution  to  a  Gaussian  distribution  in  this  case.  On  the  right  figure  we  plot  the  error  in  the  mean  and  the 
variance.  We  see  from  the  semi-log  plot  that  as  the  order  of  Hermite-chaos  expansion  increases,  the  error  of  mean 
and  variance  decreases  exponentially  fast.  This  is  due  to  the  fact  that  the  chaos  expansion  is  a  spectral  expansion 
in  the  random  space.  Similar  exponential  convergence  rate  has  been  demonstrated  for  first-order  ODE  for  various 
Askey-chaos  basis  in  [15].  It  is  worth  noting  that  if  the  appropriate  chaos  basis,  in  this  case  the  Hermite-chaos 
corresponding  to  the  Gaussian  inputs,  is  not  chosen,  the  exponential  convergence  may  not  be  maintained  [15]. 


time  p 

Figure  2:  Solution  with  Gaussian  random  inputs  by  Hermite-chaos;  Left:  Solution  of  the  dominant  random  modes, 
Right:  Error  convergence  of  the  mean  and  the  variance. 


4  Incompressible  Navier- Stokes  Equations 

In  this  section  we  present  the  solution  procedure  for  solving  the  stochastic  Navier-Stokes  equations  by  generalized 
polynomial  chaos  expansion.  The  randomness  in  the  solution  can  be  introduced  through  boundary  conditions,  initial 
conditions,  forcing,  etc.. 

4.1  Governing  Equations 

We  employ  the  incompressible  Navier-Stokes  equations 


V  •  u 


(18) 


du 

—  +  (u  •  V)u 
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0, 

-Vn  +  i2e_1V2u, 


(19) 


where  II  is  the  pressure  and  Re  the  Reynolds  number.  All  flow  quantities,  i.e.,  velocity  and  pressure  are  considered 
as  stochastic  processes.  A  random  dimension,  denoted  by  the  parameter  0,  is  introduced  in  addition  to  the  spatial- 
temporal  dimensions  (x,  t),  thus 

u  =  u(x,t;0);  II  =  II(x,  t;9).  (20) 

4.2  Chaos  Expansion 


We  apply  the  generalized  polynomial  chaos  expansion,  or  the  Askey-Chaos  (2),  to  these  quantities  and  obtain 


u(x, t;  0)  =  J2ui(x,t)$i(€(0));  n (x,£;0)  =  Eni(x’  t)$i(£(0)). 


i= 0 


2= 0 


Substituting  (21)  into  Navier-Stokes  equations  we  obtain  the  following  equations 

p 

Yy- Ufa  t)$j  =  o, 

2= 0 

E  dy t]  ^  +  E  Ef(u*  •  v)uj)]^^i  =  -  E  vn*(x’  E  y2w$i- 


(21) 


(22) 

(23) 


2=0  2=0  j— 0  2=0  2=0 

We  then  project  the  above  equations  onto  the  random  space  spanned  by  the  basis  polynomials  {4b}  by  taking  the 
inner  product  of  above  equation  with  each  basis.  By  taking  <  - ,  <4>fc  >  and  utilizing  the  orthogonality  condition  (3), 
we  obtain  the  following  set  of  equations: 

For  each  k  =  0, . . .  P, 


=  0,  (24) 

=  VII,,  -  /?<  :  ;3V2u,...  (25) 

where  e^-fc  =<  <f\;<I)J<I>A.  >.  The  set  of  equations  consists  of  (P  +  1)  system  of  deterministic  ‘Navier-Stokes-like’ 
equations  for  each  random  mode  coupled  through  the  convective  terms. 

4.3  Numerical  Discretization 

Discretization  in  space  and  time  can  be  carried  out  by  any  conventional  method.  Here  we  employ  the  spectral//ip 
element  method  in  space  in  order  to  have  better  control  of  the  numerical  error  [29] .  The  high-order  splitting  scheme 
together  with  properly  defined  consistent  pressure  boundary  conditions  are  employed  in  time  [30] .  In  particular,  the 
spatial  discretization  is  based  on  Jacobi  polynomials  on  triangles  or  quadrilaterals  in  two-dimensions,  and  tetrahedra, 
hexahedra  or  prisms  in  three-dimensions. 


V  •  ufc 

^  +  <  ^2  >  EEe^i(Ui  •  vh-)] 

2=0  j= 0 
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4.4  Post-Processing 


After  solving  the  deterministic  expansion  coefficients,  we  obtain  the  analytical  form  (in  random  space)  of  the  solution 
process.  It  is  possible  to  perform  a  number  of  analytical  operations  on  the  stochastic  solution  in  order  to  carry  out 
other  analysis  such  as  the  sensitivity  analysis.  The  mean  solution  is  contained  in  the  expansion  term  with  index  of 
zero.  The  second-moment ,  i.e. ,  the  covariance  function  is  given  by 

-Ruu(xi,  ti;  x2,  t2)  =  <  u(xi,ti)  -  u(xi,ti),u(x2,t2)  -  u(x2,t2)  > 

p 

=  ^  [m(xi,ti)ui(x2,t2)  <  >]  .  (26) 

»= 1 

Note  that  the  summation  starts  from  index  [i  =  1)  instead  of  0  to  exclude  the  mean,  and  that  the  orthogonality  of 
the  Askey-Chaos  basis  {4>i}  has  been  used  in  deriving  the  above  equation.  Similar  expressions  can  be  obtained  for 
the  pressure  field. 

Implementation  details  and  verifications  of  the  stochastic  Navier-Stokes  solver  can  be  found  in  [16]. 

5  Flow- Structure  Interactions 

In  this  section  we  consider  two-dimensional  vortex-induced  vibrations  of  an  elastically-mounted  circular  cylinder 
subject  to  stochastic  inputs.  The  computational  domain  is  shown  in  figure  3  where  the  circular  cylinder  with  unit 
diameter  ( D  =  1)  is  located  at  the  origin.  The  size  of  the  domain  is  [—15,25]  x  [—9,9].  There  are  412  triangular 
elements  and  sixth-order  Jacobi  polynomial  are  found  to  be  sufficient  to  resolve  the  flow  in  the  range  Re  <  200.  The 
Reynolds  number  is  defined  as  Re  =  UqoD /v,  where  Uoa  is  the  inflow  and  v  the  kinematic  viscosity. 


Figure  3:  Schematic  of  the  domain  for  flow  past  an  elastically  mounted  circular  cylinder. 
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5.1  Structure  Problem 


In  this  paper  we  will  focus  on  the  cross-flow  displacement  of  the  cylinder,  i.e.,  the  cylinder  is  free  to  move  in  the 
p-direction  but  not  in  the  x-direction.  For  a  linear  structure,  the  governing  equation  is  the  second-order  ordinary 
differential  equation 

£+l%+K,.m  m 

where  p,  b  and  K  are  the  mass,  damping  and  stiffness  of  the  cylinder,  and  the  natural  frequency  of  this  system  is 


Lon  =  \J I\ / p.  For  clarity,  we  re-write  equation  (27)  in  the  same  form  as  in  (11): 

_+c-  +  t,  =  /((), 


(28) 


where  c  =  b/p ,  k  =  K/p  and  /(f)  =  F(t)/p.  The  external  force  /(f)  comes  from  the  flow  and  we  incorporate 
uncertain  components  in  c  and  k  in  the  following  simulations. 

5.2  Transformed  Navier-Stokes  Equations 


To  couple  the  flow  with  moving  boundaries  of  the  structure,  one  can  employ  Arbitrary  Lagrangian-Eulerian  (ALE) 
method.  Although  general,  this  approach  is  computationally  expensive  so  we  consider  a  boundary-fitted  coordinate 
approach  for  the  specific  problem  we  solve  here.  By  attaching  the  coordinate  system  to  the  cylinder,  the  cylinder 
appears  stationary  in  time  (with  respect  to  that  coordinate  system).  Following  [31],  we  define  two  coordinate  systems: 
(x',y',t')  and  ( x,  y,t ),  where  (x',y',tr)  is  the  original  coordinate  system  and  (x,y,t)  is  the  transformed  one.  The 
mapping  between  the  two  systems  is 

{x  =  x', 

y  =  y' 

t  =  t’. 


In  two-dimensional  flow,  this  simply  reduces  to  the  v  velocities  being  shifted  by  the  reference  frame  velocity, 


v  =  v'  — 
P  =  v’  ■ 


dr] 
dt'  ’ 


It  is  worth  noting  that  this  mapping  is  stochastic  when  the  cylinder  motion  is  random  and  needs  to  be  represented 
by  the  chaos  expansion  as  well. 

The  incompressible  Navier-Stokes  equations  are  transformed  into: 


where  Ax  =  0  and  Ay  = 


d2  ij 
dt2  ’ 


V  •  u 


<9u 

—  +  (u  •  V)u 


0, 

— VII  +  f?e-1  V2u  +  A(f), 


(29) 

(30) 
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In  the  following  simulations,  we  assume  the  damping  c  and  stiffness  k  in  equation  (28)  to  be  random  variables. 
Then  the  structure  response  becomes  a  random  process,  so  does  the  psuedo-forcing  A(t)  in  the  transformed  Navier- 
Stokes  equations.  This,  in  turn,  makes  the  flow  field  random,  which  exerts  a  stochastic  dynamic  forcing  f(t)  back 
onto  the  cylinder.  The  entire  coupled  system  then  becomes  stochastic.  The  same  expansion  procedure  as  in  section 
4  is  employed,  with  the  psuedo-forcing  A (f)  and  the  mapping  expanded  appropriately,  too. 

5.3  Numerical  Results 

We  assume  that  the  Reynolds  number  is  fixed  at  Re  =  100,  and  we  also  assume  that  the  input  parameters  of 
the  cylinder  are  uncertain,  i.e.  c  =  c  +  crct;i  and  k  =  k  +  <Jk£,2,  where  and  £2  are  two  independent  standard 
Gaussian  random  variables  with  zero  mean  and  unit  variance.  The  mean  and  standard  deviation  of  c  and  k  are  set 
as  (c,  ac)  =  (0.1,0.01)  and  (k,tjk)  =  (1.0, 0.2),  respectively.  We  choose  k  =  1.0  such  that  the  natural  frequency  of 
the  oscillator  is  close  to  the  frequency  of  the  vortex  shedding  of  the  fixed  cylinder  at  Re  =  100,  and  the  cylinder 
response  is  maximized.  According  to  the  uncertain  inputs,  we  employ  the  two-dimensional  (n  =  2)  Hermite-chaos, 
the  corresponding  Askey-chaos  for  Gaussian  inputs  as  shown  in  table  1,  as  the  trial  basis  in  random  space.  A  third- 
order  Hermite-chaos  ( p  =  3)  is  employed  which  results  in  a  10-term  chaos  expansion  (P  =  9  according  to  equation 
(16)).  The  fluid  forces  on  the  cylinder  are  computed  using 

F  =  [— nil  +  Pe_1(Vu  +  VuT)  •  n]  ds, 

where  n  is  the  outward  normal  on  the  cylinder  and  ds  is  the  arc  length  on  the  surface  of  the  cylinder.  The 
corresponding  force  coefficients  are  computed  by  non-dimensionalizing  the  forces  with  the  fluid  density  pj ,  free- 
stream  velocity  Uoo  and  the  cylinder  diameter  D: 

C  =  F°  r  =  Fl 

D  \pfDUl V  L  \pfDUl' 

In  figure  4  we  plot  the  time  evolution  of  the  first  few  coefficients  of  the  dominant  random  modes  of  the  nondi- 
mensional  cross-flow  displacement  ( y/D )  and  the  lift  coefficient  (Cl),  together  with  the  deterministic  solution.  We 
see  that  due  to  the  effective  diffusion  of  the  randomness,  the  mean  response  of  y/D  has  smaller  amplitude  compared 
to  its  deterministic  counterpart.  The  first  and  second  random  modes,  as  shown  in  the  figure,  correspond  to  the 
Gaussian  part  of  the  response. 

In  figure  5  we  show  the  time  evolution  of  the  variances  of  the  cross-flow  displacement  y/D  and  lift  coefficient  Cl- 
We  see  that  the  variance  peaks  at  the  early  transition  stage  before  it  settles  to  the  asymptotic  periodic  state.  The 
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Figure  4:  Dominant  random  modes  of  the  cylinder  motion:  Upper:  Modes  of  the  cross-flow  displacement  y / D,  Lower: 
Modes  of  the  lift  coefficient  Cl- 
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peak  value  is  2  ~  3  times  larger  than  that  of  the  final  periodic  state.  This  suggests  that  the  system  responses  to  the 
uncertain  inputs  are  important  in  the  early  transition  stage  and  also  non-negligible  in  the  final  asymptotic  state. 


Figure  5:  Variance  of  the  cylinder  motion:  Upper:  Variance  of  the  cross-flow  displacement  y/D,  Lower:  Variance  of 
the  lift  coefficent  Cl- 

Figure  6  shows  the  instantaneous  contours  of  the  rms  of  the  vorticity  field  at  t  =  600  (nondimensional  time  units) 
corresponding  to  more  than  100  shedding  cycles.  The  location  of  the  cylinder  is  not  at  the  center  as  shown  in  the 
figure.  It  is  interesting  that  the  regions  with  the  largest  uncertainty  are  regions  of  the  most  importance  from  the 
fluid  dynamical  point  of  view,  i.e.  the  shear  layer  and  near-wake  but  not  the  far-fielcl. 

In  figure  7  the  instantaneous  pressure  distribution  along  the  surface  of  the  cylinder  is  shown.  Here  9  is  the  angle 
of  the  location  on  the  surface  with  9  =  0  the  rear  stagnation  point  and  9  =  ir  the  front  stagnation  point.  The 
error-bar  curve  is  centered  at  the  mean  of  the  stochastic  pressure  solution  and  the  length  of  the  bars  indicates  two 
standard  deviations  around  the  mean  (i.e.,  one  above  and  one  below  the  means).  For  comparison,  the  determinstic 
pressure  distribution  at  the  same  instance  is  shown  as  well.  The  difference  compared  with  stochastic  mean  solution 
is  noticeable;  for  the  chosen  magnitudes  of  variance  of  the  stochastic  inputs,  the  deterministic  signal  remains  inside 
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Figure  6:  Regions  of  uncertainty:  Instantaneous  rms  of  vorticity. 
the  ‘envelope’  of  the  stochastic  solution. 

6  Summary  and  Discussion 

We  have  developed  a  stochastic  spectral  method  to  model  uncertainty  and  its  propagation  in  flow  simulations.  More 
specifically,  we  have  generalized  the  original  polynomial  chaos  idea  of  Wiener  and  proposed  a  broader  framework, 
i.e.  the  Askey-Chaos,  which  includes  Wiener’s  Hermite-Chaos  as  a  subset.  Numerical  examples  were  presented  for 
relatively  simple  systems,  such  as  a  second-order  ordinary  differential  equation  and  a  more  complicated  flow-structure 
interaction  problem. 

The  method  we  developed  here  is  general  and  can  also  be  applied  to  model  uncertainty  in  the  boundary  domain, 
e.g.  a  rough  surface,  in  the  transport  coefficients,  e.g.  the  eddy  viscosity  in  large  eddy  simulations,  and  other 
problems.  It  provides  a  formal  procedure  for  constructing  a  composite  error  bar  for  CFD  applications,  as  proposed 
in  [32],  that  includes,  in  addition  to  the  discretization  errors,  contributions  due  to  imprecise  physical  inputs  to  the 
simulations. 

As  regards  efficiency,  a  single  Askey-Chaos  based  simulation,  albeit  computationally  more  expensive  than  the 
corresponding  deterministic  solver,  is  able  to  generate  the  solution  statistics  in  a  single  run.  In  contrast,  for  the 
Monte  Carlo  simulation,  tens  of  thousands  of  realizations  are  required  for  converged  statistics,  which  is  prohibitively 
expensive  for  most  CFD  problems  in  practice. 
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9/ji 

Figure  7:  Instantaneous  pressure  distribution  along  the  surface  of  the  cylinder;  Error-bar:  stochastic  solution,  Dashed 
line:  Deterministic  solution. 
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